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The radiation dynamics of a dipole antenna embedded in a Photonic Crystal are modeled by an 
initially excited harmonic oscillator coupled to a non-Markovian bath of harmonic oscillators repre- 
senting the colored electromagnetic vacuum within the crystal. Realistic coupling constants based 
on the natural modes of the Photonic Crystal, i.e., Bloch waves and their associated dispersion 
relation, are derived. For simple model systems, well-known results such as decay times and emis- 
sion spectra are reproduced. This approach enables direct incorporation of realistic band structure 
computations into studies of radiative emission from atoms and molecules within photonic crystals. 
We therefore provide a predictive and interpretative tool for experiments in both the microwave and 
optical regimes. 
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I. INTRODUCTION 

Photonic crystals (PCs) have been the subject of in- 
tensive research over the past decade [[j] . These are fab- 
ricated periodic dielectric arrays that employ a combi- 
nation of (i) Mie scattering from individual elements of 
the array, and (ii) Bragg scattering from the periodic lat- 
tice, to induce a band structure for photon propagation. 
This band structure is, in many ways, analogous to elec- 
tronic band structure in a semiconductor. Through a 
judicious selection of materials and of the periodicity of 
the lattice, the photonic dispersion relation and the as- 
sociated electromagnetic (EM) mode structure of a PC 
can be adapted to a variety of device applications. The 
most dramatic modification of the photon dispersion oc- 
curs when the linear propagation of a photon in a PC 
is prohibited in all directions for a range of frequencies, 
giving rise to a complete photonic band gap (PBG). 

The radiative dynamics of an optically active material 
placed within or near a PC can be dramatically modified 
from free space. This is a result of the "colored" electro- 
magnetic reservoir provided by the solutions to the elec- 
tromagnetic field equations within a PC. In the optical 
domain, theoretical studies of atomic transitions coupled 
to the EM modes of a PC with an optical PBG predict 
a number of novel quantum optical phenomena. These 
phenomena include: the suppression or enhancement of 
spontaneous emission and the associated fractional lo- 
calization of light near radiating atoms [^H; rapid all- 
optical switching and anomalous superradiant emis- 
sion, as well as low-threshold lasing near the edge of a 
PBG |,§. Microfabrication of PCs with complete PBGs 
at optical wavelengths has proven to be a difficult task, 
both because the lattice periodicity should be compara- 
ble to the wavelength of the light under consideration, 
and because a high dielectric contrast between the ele- 
ments of the lattice is required. Fortunately, recent ad- 
vances in microlithography Hand in semiconductor in- 



filtration in colloidal crystals |g{ have produced materi- 
als with significant pseudo-gaps in their photonic band 
structure B]. The development of materials with com- 
plete PBGs in the optical regime appears imminent. 

High-quality PBG materials at microwave frequencies 
have been available for some time |10|. Sizeable band- 
gaps with center frequencies ranging from a few GHz up 
to 2 THz have been reported; these crystals have thus 
proven the soundness of the concept of the PBG. Mi- 
crowave PBG materials may be relatively easily manu- 
factured using micro-machining techniques, and are cur- 
rently of interest for applications such as the shielding of 
human tissue from microwave radiation, and for improv- 
ing the radiation characteristics of microwave antennae. 
Although PBG materials at microwave frequencies have 
been extensively studied, the behavior of radiating dipo- 
lar antennae embedded in microwave PCs has not re- 
ceived the same degree of attention. This is despite the 
fact that such antennas would share many properties in 
common with atomic emission in a PC. In the microwave 
domain, a dipole antenna could take the form of an elec- 
trically excited metallic pin with a high Q (quality) fac- 
tor. 

The radiative dynamics of the above system can be de- 
scribed by a charged, one-dimensional simple harmonic 
oscillator (SHO). Such an electric dipole oscillator can 
also provide an excellent description of the radiation of 
single or multiple two-level atoms in the optical domain. 
This description is valid provided that the total excita- 
tion energy of the atoms is well below an energy where 
saturation (nonlinear) effects become important. More- 
over, the radiation reservoir can itself be modeled as a 
bath of many independent SHOs: Radiative damping 
arises from a linear coupling between the system SHO 
and the large number of reservoir oscillator modes. The 
similarities between the microwave and optical systems, 
coupled with the mature state of microwave technology, 
suggests that many of the predicted effects for atomic 
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dipoles in the optical domain could be realized and stud- 
ied first in the microwave domain. 

Analytical techniques exist for treating certain forms 
of coupling between the dipole and reservoir for cer- 
tain modal distributions of the reservoir. However, 
PCs present coupling distributions and spectral prop- 
erties which defy analytical methods. This is due to 
the presence of a restricted and rapidly-varying reservoir 
mode distribution, which renders invalid the usual Born- 
Markov type approximation schemes for the system- 
reservoir interaction. To obtain accurate results, we solve 
the system numerically for a large, but finite, number 
of oscillators in the reservoir by discretizing the modes 
of the reservoir following the approach of Ullersma . 
In dealing with our system, there are crucial issues con- 
cerning obtaining the correct coupling strength between 
the oscillator and the reservoir modes, as well as in em- 
ploying the proper renormalization and mode sampling 
in numerical simulations. When these criteria are satis- 
fied, the SHO method comprises a powerful approach to 
treating radiative dynamics. 

Here, we develop a rigorous quantitative treatment of 
the radiative dynamics of an electric dipole oscillator cou- 
pled to the electromagnetic reservoir within a model PC. 
In the process, we provide a sound theoretical basis for 
this and other approaches p2| to non-Markovian radia- 
tive dynamics which involve the discretization of a model 
electromagnetic reservoir. Additionally, we show how our 
method can be applied to realistic PC's with complicated 
dispersion relations and EM mode structures. The pa- 
per is organized as follows. In Section ||, we develop a 
classical field theory for electromagnetic field modes in 
PCs, and we derive the coupling constants for the cou- 
pling between a radiating dipole and these Bloch modes. 
This leads to the Hamiltonian of the coupled system and 
the associated equations of motion. Renormalization is- 
sues arising from the non-relativistic nature of our the- 
ory are discussed in Section III, whereas Section IV de- 
scribes the discretization of the reservoir and the numer- 
ical solution of the equations of motion. In Section V, 
these techniques are applied to a highly computationally 
challenging model, that of a three-dimensional, isotropic 
dispersion relation with a complete PBG. The demon- 
stration of fractional localization and related phenomena 
validates the SHO approach to modeling radiative dy- 
namics in PCs. In Section VI we summarize the results 
and emphasize the possibilities for testing these predic- 
tions experimentally in the microwave domain. The two 
appendices are concerned with the details of the field the- 
ory for the PC and with the details of the model of the 
one-sided, isotropic PBG, respectively. 

II. CLASSICAL FIELD THEORY 

In this section, we derive the equations governing the 
dynamics of a radiating dipole oscillator located inside a 



PC. Typically the equation of motion for a damped os- 
cillator, with time-dependent coordinate q(t), is written 
as the second-order differential equation 

«(*)+7«(*)+<4z(t) = F(t). (1) 

Here, we have introduced a damping constant 7, the 
natural frequency uq and the driving field F(t) for the 
amplitude q of the linear oscillator. For instance, for a 
freely oscillating RLC circuit with ohmic resistance R, 
capacitance C and inductance L, we have 7 = R/L, 
cj% = 1/LC, F(t) = 0, and q(t) is the electric charge. Eq. 
(0) is, however, not the most general way of incorporat- 
ing damping into the equations of motion for a harmonic 
oscillator. This description can break down if, for exam- 
ple, there is a suppression of modes in the reservoir to 
which the dipole oscillator can couple. Such a suppres- 
sion of modes is a feature of the EM reservoir present 
in a PC. A more general description of damping forces 
acting on the harmonic oscillator therefore requires a pre- 
cise knowledge of the mode structure of its environment, 
and the corresponding coupling of the system oscillator 
to these modes. In the case of a radiating dipole located 
in a PC, it is then appropriate to model its emission 
dynamics with a SHO coupled to a reservoir of SHOs. 
The essential difference between the vacuum and a PC 
is then contained in the spectral distribution, or density 
of states (DOS), of the reservoir oscillators, and in the 
coupling constants between the reservoir modes and the 
system oscillator. 

The characterization of the reservoir is carried out in 
detail in Appendix A; here we only summarize the salient 
results. Given a radiating dipole with a natural frequency 
u>o, we obtain the classical Hamiltonian 

H = i?dip + -Hires + H ct + Hi n t ■ (2) 

The first term on the right-hand side of the Hamiltonian 
is the energy of the dipole oscillator itself, 

#di P = ^oH 2 . (3) 

The natural frequency of the isolated oscillator is ujq, 
and £ is a constant with the dimension of energy x time. 
This permits us to write the energy of a SHO in units 
of its natural frequency w, i.e., E(u>) = ^lu. The system 
oscillator's complex amplitude is given by the dimension- 
less, time-dependent quantity a, defined with respect to 
the coordinate q(t) of Eq. (Q) as 

aw-^w+TS^)) • (4) 

The next term in the Hamiltonian (0) corresponds to 
the free evolution of the radiation reservoir, which is 
modeled as a bath of independent SHOs, 

H tea = J2t^\P»\ 2 - (5) 



2 



The natural electromagnetic modes of the PC are Bloch 
modes (see Appendix A), labeled with the index /i = 
(nk), where n stands for the band index and k is a re- 
ciprocal lattice vector that lies in the first Brillouin zone 
(BZ). Their dispersion relation, is different from the 
vacuum case, and may have complete gaps and/or the 
corresponding density of states may exhibit apprecia- 
ble pseudogap structure, the manifestation of multiple 
(Bragg) scattering effects in periodic media. 

As we are working within the framework of a non- 
relativistic field theory, we have introduced a mass renor- 
malization counter term H ct — — fA|a| 2 that cancels un- 
physical UV-divergent terms jl3 14|. The quantity A is 
specified in Section III. 

The interaction between the oscillator and the reservoir 
is given by a linear coupling term. As the oscillator fre- 
quency is quite large, and the effective linewidth of the 
oscillation is relatively small, it is possible to simplify 
the interaction by applying the rotating-wave approxi- 
mation. In this approximation, couplings in the Hamil- 
tonian of the form a* (3* and its complex conjugate are 
neglected, as these terms oscillate very rapidly compared 
to the terms of the type a*/3 M and its conjugate. Hence, 
the interaction Hamiltonian can be expressed as 



Hi n t = ( a * 9*p, P» -"S^J 



(6) 



In the case of a point dipole, i.e., when its spatial ex- 
tent a is much smaller than the wavelength corresponding 
to its natural frequency, Ao = 2ttu)q/c, the coupling con- 
stants can be derived from (i) the magnitude of the 
dipole moment, d(t) = aq(t), located at ro, and (ii) the 
dipole orientation, d, relative to that of the Bloch modes, 



9p, = 9n(ro) = ac 



Lujoujfj, 



(7) 



This dependence of the coupling constant on the dipole's 
precise location within the PC is the second essential 
difference from the free-space case. As shown in Refs. 
JDiLfni, this position dependence may be quite strong, 
thus making its incorporation a sine qua non for any 
quantitative theory of of radiating antennae or fluores- 
cence phenomena in realistic PCs. 

The emission dynamics can be evaluated from the Pois- 
son brackets of the oscillator amplitudes and their initial 
values, a(0) = 1 and /3 P (0) = (V/i). Our choice of a(0) 
and P/i(0) corresponds to the initial condition of an ex- 
cited dipole antenna and a completely de-excited bath. 
The only non-zero Poisson brackets are 



{a,a*} = 



(8) 



Eqs. (|2|), (0) and (||), together with the initial values 
for the oscillator amplitudes, completely determine the 
emission dynamics of a radiating dipole embedded in a 



PC. In the following sections, we solve the correspond- 
ing equations of motion. This task is complicated by 
the nature of the reservoir's excitation spectrum: as dis- 
cussed, the non-smooth density of states prohibits the use 
of a Markovian approximation and its appealing simpli- 
fying features Instead, we have to revert to a 
solution of the full non-Markovian problem. This is ac- 
complished by firstly rearranging the reservoir modes in 
a manner more suitable to both analytical as well as nu- 
merical solutions, and subsequently solving the equations 
of motion. In what follows, we bridge the gap between 
previous studies of simplified model dispersion relations 
and band structure computations p5|Jl7| |. 
Although we will formally develop our theory for an 
LC circuit in a microwave PC, we emphasize that the 
formalism applies equally well to a semiclassical Lorentz 
oscillator model of an excited two- level atom, i.e., an 
electron with charge e and mass m which is bound to 
a stationary nucleus, for which the energy of excitation 
is well below that required for saturation effects to be- 
come relevant. The oscillator coordinate q(t) may then 
be identified with the deviation of the electron's position 
from its equilibrium value, 7 is the inverse life time of 
the excited state, and uiq denotes the frequency for tran- 
sitions between excited and ground state of the two-level 
atom. This corresponds to making the substitutions: 



L^m, (Lq) p, £ -> h, 
where h = 2irh is Planck's constant. 



(9) 



III. PROJECTED LOCAL DENSITY OF STATES, 
MASS RENORMALIZATION AND LAMB SHIFT 

From the Hamiltonian (j^) we derive the equations of 
motion for the amplitudes 

i(t) = -t (uo — A) (10) 



a 



0n(t) = -iLo^fi^it) +5 M a(t), 



(11) 



for which we seek a solution with initial conditions a(0) = 
1 and /3^(0) = (V/i). Our formalism however requires 
that we first determine the mass renormalization counter 
term A. This is most conveniently done in a rotating 
frame with slowly varying amplitudes a(t) and b(t), de- 



fined as a(t) 
tively: 



a(t)e- iuJot and (3(t) = b{t)e 



respec- 



o(*) = -%t E 9l e^o-^X 6 P (t) + »Ao(t) 



-z(cj, 



(12) 



b(t)=g„e-^ u °-"»> t a(t). (13) 

Conversely, Eqs. ( |l2j ) and ( |l3| ) comprise a stiff set of dif- 
ferential equations making their solution a difficult task. 
Numerical solution of the problem is more easily per- 
formed in the non-rotating frame, to which we return in 
Sect. IV. 
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Eq. ( |l3| ) may be formally integrated, 

b lt (t)=g ll ( dt' e-^-"^' a{t'), (14) 
Jo 

and inserted into Eq. (|l^) to yield 

poo 

a(t) = - dt'G(t-t')a(t') + iAa(t), (15) 



where the Green function G(t) contains all the informa- 
tion about the reservoir and is the subject of our studies 
for the remainder of this section. It is defined as 



(16) 



Here, O(r) denotes the Heaviside step function, which 
ensures the causality of G(r). We now proceed to evalu- 
ate G(t) for the form of the coupling constants <? M given 
in Eq. (|7[). To this end, we introduce the projected local 
DOS (PLDOS) N(r ,d,co) through 



N(f ,d,cj)=J2 



d 3 k 



8{u-u n% )\d-E n% {f )\\ 

BZ \ Z7r ) 



(17) 



where we have replaced the symbolic sum over /i by its 
proper representation as a sum over bands plus a wave 
vector integral over the BZ. With these changes, we may 
rewrite G(r) compactly as 



G(t)=PQ(t) / dco 
Jo 



N(f ,d,u;) e4(wo _ w)r 



(18) 



Here, we have abbreviated /3 = (na 2 c 2 )/ (Lloq). Eq. (|l§| ) 
makes more explicit what we have argued before: The 
spontaneous emission dynamics of active media in Pho- 
tonic Crystals are completely determined by the PLDOS, 
N(rQ, d, lo). As the PLDOS may be drastically different 
from location to location within the Wigner-Seitz cell of 
the PC (l^Jl^], it is imperative to have detailed knowl- 
edge about where in the PC the dipole is situated in order 
to understand and predict the outcome of corresponding 
experiments. 

One additional point deserves special attention: the 
total DOS, N(u), is related to the local DOS via 

N(w) = ± J d 3 r J dn d e p (r)N(r,iio) 



J v J MiNftlw), 

where V is the volume of the Wigner-Seitz cell, and 
J dilj is the average over all possible orientations of the 
dipole. Strictly speaking, it is not possible to base con- 
clusions about the radiation dynamics on the total DOS. 
This is a direct consequence of the fact that the natural 



modes of PCs are Bloch waves rather than plane waves as 
in free space. Depending on the band index, these Bloch 
modes prefer to "reside" predominantly in either low or 
high dielectric index regions (so-called air and dielectric 
bands respectively) . Only in the case of very low index 
contrast ("nearly free photons") may the total DOS be 
viewed as a reliable guide to interpreting radiative dy- 
namics within a PC. The total DOS is, nevertheless, an 
adequate rule-of-thumb estimator. 

From Eq. (0) we can now obtain the Fourier transform 
of the Green function, G(f2 — loq), centered around the 
atom's bare transition frequency too: 

/>oo 

G(fi-wo) = / dte< Q -" ^G{t) 



Q — lo 



where p stands for the principal value. 

For large lo, we have N(fo, d,Lo) cx lo 2 . The imaginary 
part of G(Q — luq) apparently contains a linear divergence 
in the UV. This divergence is to be expected for a non- 
relativistic theory, analogous to the problem of sponta- 
neous emission in vacuum |l5|], and is removed from the 
theory by using the mass counter renormalization term 
A, as first pointed out by Bethe [Q. Consequently, we 
decompose the imaginary part of G(f2 — luq) into 



3 (G(f! - Wo )) ~ -(A + <5 vac + 5 a ) , 
where we have used the notation: 



(19) 



A = (3 duo 
Jo 

vac — r> o 

S = 

7T 2 C 3 jq 



duo p 



duo p 



1 



LUq — LO 
1 



k LUq — LO / 

N(r Q ,d,Lo) ~ N^ito) 

X 2 ' 

uo z 

Here, we have performed a Wigner-Weisskopf-type ap- 
proximation on the vacuum and anomalous Lamb shifts 
0, <5vac ancl <^a> respectively. This approximation is jus- 
tified by the fact that, despite its highly non-Mar kovian 
nature, a radiating dipole in a PC is still a weak coupling 
problem, as can be seen, for instance, by estimating the 
coupling constant 



2tt 



(20) 



is the 



in the Lorentz oscillator model. Here, V ~ 
volume of the Wigner-Seitz cell of the PC (a is the cor- 
responding lattice constant) and do = eao is the oscil- 
lator's dipole moment for the elementary charge e and 
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Bohr atomic radius, a,Q. At optical frequencies (ui « 10 15 
s _1 ), a silicon inverted opal has a PBG at the fre- 
quency aui/2irc « 0.8, so that we obtain 10 -7 < g/uuo < 
1CP 6 <C 1, thus justifying our Wigner-Weisskopf approx- 
imation. As a consequence, we must treat the real part 
of G(f2 — ujo) exactly, but are still allowed to tackle the 
imaginary part of G(Q — luq) using standard perturba- 
tion methods of QED. In addition, we have introduced 
the vacuum or free-space DOS A( vac )(w) = uj 2 /(tt 2 c 3 ), 
and a cutoff frequency Q c 3> luq, which is chosen large 
enough that the results of the following analysis remain 
independent of the precise value of f2 c . In a Lorentz os- 
cillator model, for instance, O c can be identified with the 
Compton frequency fl c ~ mc 2 /%, as u> > uj c probes the 
relativistic aspects of the oscillating charge, which are 
beyond the scope of the model. 

With the foregoing analysis, we have determined the 
mass renormalization counter term A. In addition, we 
have derived an explicit expression for the anomalous 
Lamb shift S a Q which originates in the "reshuffling" 
of the reservoir's spectral weight by the PC. 



IV. DISCRETIZATION OF THE RESERVOIR 

To solve the equation of motion for the amplitude of 
the system oscillator, let us rewrite Eq. (|l5|) in a more 
explicit form: 

/■OO f-t 

d(t) = - / dujN(r ,d,uu)g 2 {cu) / dtV^-^^'Vt') 
Jo Jo 

-HAa(t), (21) 

where g 2 (co) = f3/u>, and the mass renormalization 
counter term A is given by 



A = f3 / dio 
>o 



N{r ,d,oj) 



(22) 



We remind the reader that a(0) = 1. 

We are now in a position to comment on the origin of 
the linear damping term ^q(t) that appears in Eq. (Q): 
If we consider the long time limit, i. e., t 1/loq, and 
assume that the PLDOS N(fo, d, uj) is a smooth function 
for frequencies around ujq, we can approximate the fre- 



quency integral in Eq. (^IJ) by 
t'), which leads to 



a(t) = -ja(t), 
where the decay constant is defined as 
7 = 2irf3N(r , d, uj )/oj . 



2<irpN(?a,d,uo)/ua\ S(t- 
(23) 



(24) 



This approximation is is valid only for long times relative 
to 1/oj , an d f° r a sufficiently smooth density of states. 
However, in the case of a PC, the PLDOS may have 



sharp discontinuities and gaps, thus requiring that the 
full equations of motion instead. 

To solve the integro-differential equation ( |2l| ) in a PC, 
we appeal to the literal meaning of the PLDOS as a den- 
sity of states: A(r*o,<i, oj) may be interpreted as an un- 
normalized probability density of finding a reservoir oscil- 
lator with frequency u) at position fn and orientation d. 
Consequently, we transform Eq. (|2l]) back to a system 
of coupled differential equations by employing a Monte 
Carlo integration scheme for an arbitrary function f(uj) 
according to 

du;N(r ,d,u;)f(w)~ / dio N(f , d, uj) J(uj) 
o Jo 

M 



No 
M 



where the normalization constant 



doj N(ro, d, uj) 



(25) 



(26) 



depends on the cutoff frequency, il c . There are M ^> 1 
bath oscillators, contained within a set of frequencies 
{uJi, 1 < i < M}, the frequencies of which are obtained 
by randomly sampling the interval [0, tt c ] according to 
the probability density p(r ,d,uj) — N(r Q ,d,uj)/No. 
Note that the u>i may be degenerate, as prescribed by 
p(r Q ,d,uj). 

Applying this Monte Carlo scheme to Eq. (|2l]) and 
transforming back to a non-rotating frame in order to 
avoid having to solve a numerically stiff problem, we ob- 
tain 

JV 

a{t) = -i (too - A) a(t)- i c£ftAW (27) 



$i(t) = -iuj,f3 l (t)+g l a(t), 



(28) 



where <?, = g(uji), 1 < i < M, and the mass renormaliza- 
tion counter term is evaluated up to the cutoff frequency 

tt c , i.e., A = J ° c da; Niroid, u>)/uj 2 . 

When comparing Eqs. @ an d (EH to our initial equa- 
tions of motion, Eqs. ( ^0| ) and (|ll), we observe that the 
considerations in the previous section have allowed us 
to rearrange the three-dimensional wave vector sum over 
the modes \x = (nk) into a simple one-dimensional sum 
over a set of frequencies {uii\ with a probability distri- 
bution p(rb, d, uj) that is easily determined through stan- 
dard photonic band structure computation ]l5| . In the 
following section, we give the solutions of (|27|) and ( psj ) 
for a model system which has previously been treated by 
other methods. In particular, we will demonstrate that 
known results for the radiative dynamics can be recap- 
tured and do not depend on the the value of the cutoff fre- 
quency f2 c and the number M of reservoir oscillators once 
these quantities are large enough such that all the rele- 
vant features of N(f , d, lu), are adequately represented. 
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V. NUMERICAL RESULTS FOR A MODEL 
SYSTEM 

In order to establish the validity of our approach, we 
now solve Eqs. ( p7| ) and (28) for a generic model of a 
PBG, the three-dimensional isotropic, one-sided PBG ||. 
In Appendix B, we outline the construction of the model's 
dispersion relation and how to obtain the corresponding 
model DOS, N m (ui). We note that we do not appeal to 
an effective mass approximation in the dispersion relation 
||, as is done in most treatments of band-edge dynamics. 
This allows us to recover the correct form of the large 
frequency behavior of the photon density of states. 

In Fig. 1, we show the behavior of N m (u>) as a func- 
tion of frequency for values of the relevant parameters, 
the gap size parameter r\ = 0.8 and the normalized cen- 
ter frequency Lu c a/2irc — 0.5 (sec Appendix B). The 
DOS exhibits a square-root singularity at the band edge 
bj u a/2nc = 0.6, as well as a UV divergence, N m (u>) oc co 2 , 
as oj — > oo; these are the characteristic features of this 
model. Due to the simultaneous presence of both diver- 
gences, this model clearly represents a severe numerical 
test of our approach. In order to test the method, we 
thus replace the PLDOS entering Eqs. ( ]27| ) and ( p8| ) by 
N m {u). 

In Fig. 2, we present the results of our numerical so- 
lution for the radiation dynamics of a dipole oscillator 
with frequency ujq that is coupled to the modes of a PC, 
as described by Eqs. ( p7| ) and (p8|), for various values 
of the bare oscillator frequency, woa/27rc, relative to the 
bandedge at Lo u a/2irc =0.6. The coupling strength has 
been chosen such that g(uJo) = 10~ 4 , corresponding to 
[3=10- h xluI 

Clearly visible are normal mode oscillations, also re- 
ferred to as vacuum Rabi oscillations, and the fractional 
localization of the oscillator's energy at long times near 
the photonic band-edge [||. As expected, for frequen- 
cies deep in the photonic band-gap {ujQa/2irc — 0.58), 
where the system oscillator is effectively decoupled from 
the bath oscillators, we find no noticeable decay of the 
oscillator amplitude. Deep in the photonic conduction 
band (<jJqcl/2-kc = 0.62), the system oscillator is coupled 
to a bath with a smooth and slowly-varying mode den- 
sity, as in free space. We therefore observe exponential 
decay of the oscillator amplitude, though with a time 
scale that differs significantly from that in free space. 
Due to the large value of the DOS close to the pho- 
tonic band edge, the initial decay is faster for bare os- 
cillator frequencies close to this edge than for frequencies 
deep inside the allowed photonic band. These results 
were obtained for a smooth exponential cutoff for the 
DOS around fl c a/2nc = 3.0 and M = 2.5 x 10 5 oscil- 
lators representing the modes of the PC. We also per- 
formed numerical simulations between all combinations 
of Q, c and M with values Q, c a/2irc — 3.0,6.0,9.0 and 
M = 2.5 x 10 5 ,5 x 10 5 , 10 6 and found that the numeri- 
cal values differ by at most 0.2% of the values shown in 



Fig 1. This demonstrates that, despite the presence of 
the singularities in the DOS, our approach still provides 
accurate and convergent results. 



VI. DISCUSSION 

In summary, we have developed a realistic field theory 
for an oscillating electric dipole located in a PC. The 
theory is based on the natural modes of the PC, the 
Bloch waves, and allows the direct incorporation of real- 
istic band structure calculations in order to obtain quan- 
titative results for the radiation dynamics of the dipole 
antenna. We have shown how the theory must be renor- 
malized in order to account for unphysical divergences 
and have identified the classical analogue of the Lamb 
shift of the dipole's natural radiation frequency. Finally, 
we have developed a reliable numerical scheme based on 
a probability interpretation of the PLDOS that solves the 
equations of motion for the dipole oscillator coupled to 
the electromagnetic mode reservoir of the PC. 

The viability of this approach was demonstrated for 
an isotropic model DOS for which we have derived well- 
known results for radiating atomic systems |^] in the con- 
text of a radiating classical dipole. The model considered 
contains two divergences, one square-root-divergence at 
the photonic band edge and a quadratic UV-divergence, 
and therefore clearly comprises the most serious test 
of our approach. More realistic models of a three di- 
mensional photonic band-edge take into account the 
anisotropy of the BZ, and therefore do not suffer from 
a band-edge singularity ||. As a result, our formalism is 
clearly more than capable of treating more realistic de- 
scriptions of the electromagnetic reservoir within a PC. 

Though we have developed our theory for an LC cir- 
cuit in a microwave PC, we have pointed out in Section 
II that the formalism applies equally well to a semiclassi- 
cal Lorentz oscillator model of an excited two-level atom. 
Therefore, our approach is applicable to both microwave 
antennas and to optical atomic transitions. However, 
technological constraints suggest that microwave experi- 
ments will likely be easier to perform than optical exper- 
iments involving single atoms. As discussed, an appro- 
priate microwave antenna could, for example, take the 
form of a high-Q metallic pin placed in or near a PC. 
The pin can then be excited by a focused ultrashort laser 
pulse that generates free carriers at one end; these carri- 
ers then undergo several oscillations across the pin before 
re-establishing charge equilibrium. The resulting signal 
could be easily detected and compared with the emis- 
sion from such an antenna positioned in free space, or 
within a homogeneous sample of the dielectric material 
that makes up the backbone of the PC under considera- 
tion. 

In its own right, such a microwave system could have 
considerable applications in radio science and microwave 
technology. For example, the PBG can be used as a fre- 
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quency filter, and can be used to fine tune the bandwidth 
of a dipole emitter with a resonant frequency near the 
edge of the gap. It may also be possible to actively modify 
the photonic band structure, effectively changing the ra- 
diation pattern of a dipole emitter. A feasible scheme for 
active band structure modification has recently been pro- 
posed in the context of optical PCs |l8|, in which the PC 
is infiltrated with a liquid crystalline material whose ne- 
matic director is aligned using applied electric fields. By 
rotating the director, it was found that the band struc- 
ture could be significantly modified, and that PBGs may 
be opened and closed altogether. Similar methods may 
be applied to the case of microwave PCs. 

Although we have concentrated specifically on the lin- 
ear model, the method of coupled oscillators can be ex- 
tended to treat a nonlinear antennae , or a collection 
of two-level atoms in a regime where saturation effects 
arise. As we have shown here, this method of coupled 
classical oscillators reproduces effects normally associ- 
ated with quantum optical calculations. We expect that 
a nonlinear oscillator model will reproduce some of the 
effects studied for a single two-level atom coupled to the 
modes of a PC without the need for quantizing the field. 
However, a classical treatment would need to be aban- 
doned if multiphoton excitations are non-negligible |L2] . 
Given that multiphoton effects are difficult to observe in 
the microwave domain and even more challenging in 
the optical domain [^o| , it is reasonable to expect that a 
classical model of radiative dynamics in a PC should be 
sufficient for foreseeable experiments. 
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APPENDIX A: CLASSICAL FIELD THEORY 
FOR PHOTONIC CRYSTALS 

In this Appendix, we present a self-contained formula- 
tion of a classical field theory for the Bloch modes of a 
PC, and we develop the Hamiltonian describing the cou- 
pling of a radiating dipole couples to these modes. As 
a first step, we review the computation of dispersion re- 
lations, and of electric and magnetic field modes from 



band structure calculations E^j. We then demonstrate 
how the results of such band structure calculations can 
be used to construct the corresponding vector potentials 
and free field Hamiltonian. Finally, we derive the full 
minimal coupling Hamiltonian for a classical radiating 
dipole embedded in a PC. This may be compared to the 
formulation of a general, quantized field theory for EM 
modes in nonuniform dielectric media in terms of a nor- 
mal mode expansion in the context of quantum optics 
& 



1. Review of band structure calculations 

We develop our theory in terms of the magnetic field H 
rather than in terms of the electric or displacement fields 
because (i) V • H = and, (ii) the transverse and longi- 
tudinal components of the magnetic field are continuous 
across the dielectric boundaries. This leads to more rapid 
convergence of the relevant Fourier series expansions. 

In a three-dimensional PC, we can write the eigenvalue 
equation for the magnetic field H as 



V x k(r)Vx H) + —H = 



(Al) 



with rj p {r) the inverse of the periodic dielectric permit- 
tivity, 



e p (f) = e b + (e Q - e b ) V] S(f- n ■ A). 



E 



(A2) 



The medium is assumed to consist of a background ma- 
terial with bulk permittivity e b and a set of scatterers, 
with bulk permittivity e a . The shape of the scatterers is 
described by the function S, i. e.„ S(r) = 1 if f lies inside 
the scatterer and zero elsewhere, distributed periodically 
at positions 



(A3) 



The notation of Eq. (A2) is obtained by defining the 
matrix A = (ai S3) and Z 3 = Z <£> Z (g) Z . The di- 
electric permittivity is spatially periodic modula fi ■ A. 
The assumption of a scalar permittivity is reasonable 
for bulk materials which are not birefringent but in no 
way restricts the considerations below. Chromatic dis- 
persion effects are considered to be negligible, thus al- 
lowing the time-dependence of the permittivity to be ig- 
nored. Let us define the dual matrix B = 2n(A~ 1 ) T . 
For B — {pi 62^3), this definition leads to the orthogo- 
nality relation 



Qi ■ bj = 2tt 8ij . 



(A4) 



Whereas the points ft- A are the real space lattice vectors, 
the points m ■ B, for m 6 Z 3 are the reciprocal lattice 



7 



vectors. The inverse permittivity can be expanded in the 
dual basis as 

^irh-B-r 



V P (r) = r ?™ e " 

rheZ 3 



(A5) 



The differential equation ( |A1| ) has periodic coefficients. 
By the Bloch-Floquet theorem we can expand the mag- 
netic field as 



where u? is spatially periodic modulo A; that is, 
= uAr + n-A). 



(A6) 



(A7) 



The set {k} labeling the solutions can be restricted to lie 
within in the irreducible part of the first Brillouin zone 
(BZ), since any value of k can then be obtained through 
a combination of group transformations with respect to 
an operation from the point group of the crystal and 
translations with respect to a reciprocal lattice vector. 
We can therefore express each wavevector k as 



k 



J. ,m 



k*-T + m-B 



(A8) 



where fc» is an element of the irreducible part of the 1. 
BZ and T an element of the crystal's point group. 

Applying the Bloch-Floquet theorem, Eq. flAq), the 
magnetic field can be expanded as 



4~ = ^EE^ A 4 A ^' • (A9) 

m A— 1 

Here A is the index of polarization and the vectors 



I -m -fc,2 k + rh- B 



\k + m- B\ 



(A10) 



form an orthonormal right-handed triad. This expansion 
inserted into Eq (Al) yields an infinite eigenvalue prob- 
lem which is then solved numerically by a suitable trun- 
cation. Typically the cardinality of the set {m} is on the 
order of 10 3 [El|. For any given fc* we obtain a discrete 
set of eigenfrequencies co n ^ and corresponding cigcnfunc- 
tions H t which we label by the band index n € Af. It 
is important to note that the expression for the electric 
field can be recovered from the magnetic field via 



E t(r) = -i 

nk V / 



V x H ,-(r) 

nk V / 



(All) 



In addition, the Bloch waves obey the following orthog- 
onality relations: 



d 3 rH*Af)-Hr,(r)^S nm S(k-k') , (A12) 



d 3 re p (r)E* n% (f)-E mP (f)^S nm S(k-k') , (A13) 



where the integration is over all space in both cases. We 
are free to choose the constants of proportionality in the 
above relations, and do so in the next subsection. 



2. Free field Hamiltonian 

Based on the above considerations, we are now in a 
position to derive the general expressions for the scalar 
and vector potential, <j){r, t) and A(r, t) respectively, for 
the classical Hamiltonian of the free radiation field. We 
find that the expressions become particularly transparent 
in the Dzyaloshinsky gauge, i.e., when <f>(r, t) = 0. Then, 



E(f,t) = - 



1 dA(f, t) 



c dt 
H(f,t) = V x A(r,t) 



(A14) 
(A15) 



and the gauge condition V • [e p (r)A(r, t)j = 0, reveals 

that in a PC the natural modes of the radiation field are 
no longer transverse. This is of importance wh en quan - 
t izing the field theory Given Eqs. ([Al]), ([All]), 
(A14) and ( A15 ), it is now straightforward to derive the 
following expansion of the vector potential A(f, t) 



A{r,t) 



d 3 k /2<c 2 



[JZ O) 3 V ^nk 



(/W*) + A* r Jr)) , (A16) 



where the time evolution of the free field is described by 

= Pnk(°) e ~ lu}n%t - The field modes obe y 



nk 



(A17) 



which is the same equation as that for the electric field 
modes E n ^(r) of Eq. (All). We now choose the normal- 



ization of A r such that 

nk 



d 3 r e p (r) A n% {r) ■ A m% , (r) = S nm 5(k - k') , (A18) 



d A r 



(Vxl £ (r)).(vxl m£ ,(r) 
-^S nm S(k-k') 



(A19) 



This also fixes the normalization in Eqs. (|Al~2[) and (]A1^) . 
As a consequence, the total electric and magnetic field are 
now given by 



-..^ , v-^ f d 3 k /27t£c 2 



„ Jbz ( 27T 



d 3 k /2<c 2 



bz (2tt) 3 Y u nt 
(WO^W+^W^sW) - (A21) 
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where we have re-introduced the electric and magnetic 
field modes, E n% {f) = {uJ ni :/c) A ni {r) and H n% {r) = V x 
A ni {r), respectively. Eqs. flA20D and ( |A2lD finally lead 
us to the free field Hamiltonian 



E 



BZ 



(A22) 



The only nonzero Poisson brackets are |/3„g) = %jw. 



3. Radiating dipole embedded in a Photonic Crystal 

We consider the insertion of a point dipole into a PBG 
structure at a prescribed location rb. The free dipole 
oscillator is described by the Hamiltonian H& lv 

Hdip = ^2L +^oV=^oM 2 , (A23) 

where the dipole's natural frequency is luq = l/LC and 
the complex oscillator amplitude a is given in terms of 
the charge q and "current" Lq as a(t) = q(t) \J Lluq/2w + 
i(Lq(t)) I \/2t;Lu>o, with Poisson brackets {a, a*} = 
The point dipole couples to the electric field via its dipole 
moment d(t) — aq{t) with orientation d, which yields the 
interaction energy 



H- m t — — 



aq(t) (d-E{r ,t) 



(A24) 



In the rotating wave approximation to the interaction 
term, the final minimal coupling Hamiltonian for a radi- 
ating dipole in a PC is 

H = H dlp + H ICS + H ct + H int ■ (A25) 

Collecting all the above results we obtain 

H = ^ Q \a\ 2 + Y,^\^\ 2 + H ct 

-<Y,( a *ti^~ a ^^)- ( A26 ) 



Here, we have introduced the symbolic index fi = (nk) 
and the coupling constants 



g» =gn(ro) = ac 



LujoLlIu 



(A27) 



In addition, in Eq. ( A25 ) we have introduced a mass 
renormalization counter term, H ct — — £A \a\ 2 in or- 
der to cancel unphysical UV-divergent terms of our non- 
relativistic theory, as discussed in the main text. For 
completeness, we list here only the nonzero Poisson 
brackets and initial values for an initially excited radi- 
ating dipole coupled to the Bloch waves of a PC. This, 
together with the Hamilton function H in Eq. ( A25 ) com- 
pletely defines our problem: 



{a,a*}= {{3^,(3;} 



(A28) 



APPENDIX B: MODEL DISPERSION RELATION 
AND DENSITY OF STATES 

A particularly stringent test of our approach's ability 
to describe the dynamics of a radiating dipole in a PC 
comes from its application to a dipole coupled to a 3D 
isotropic photon dispersion model for the electromagnetic 
reservoir. In this model, the coherent scattering condi- 
tion that characterizes the photonic band edge is assumed 
to occur at the same frequency for all directions of propa- 
gation. Clearly this is not the case in a real crystal, whose 
Brillouin zone cannot have full rotational symmetry. As 
a result, the isotropic model overestimates the electro- 
magnetic modes available at a band-edge, so that, for 
example, near the upper photonic band edge at frequency 
u) u , the corresponding DOS exhibits a divergence of the 
form N(u>) cx 1 / y/uj — ui u . Conversely, for large frequen- 
cies (oj 3> io u ) the DOS will exhibit a UV-divergence, i.e., 
N(ui) oc uj 2 , as is the case in free space. More realistic 
LDOS' coming from full photonic band structure compu- 
tations do not suffer from the pathological band edge di- 
vergence of the isotropic model. However, by solving the 
model of a 3D isotropic photonic band gap, we make con- 
tact with previous results based on the isotropic model 
in the effective mass approximation || . 

Consider a ID photonic dispersion relation in the ex- 
tended zone scheme. In order to describe a PBG at wave 
number k Q with central frequency oj c = cfc = (lu u + lui)/2 
and upper and lower band edge at uj u and lui , respectively, 
we use the following Ansatz 



, oj+ + c+J(k- k ) 2 +^ 2 forfc>0 

uo(k) = { + + V '± . (Bl) 

w_ + c_ J (k — ko) 2 + 7_ for k < 



Using the requirements Lo(k — 0) = 0, ui(k = ko — 0+) = 
u>i, ui(k = fco+0+) = lo u , dkuj{k — 0) = dkLo(k — > oo) = c, 
and <9fcO>(fc = k — + ) = <9fcw(fc = k + 0+) = 0, the 



unknown parameters in (Bl) can easily be expressed in 
terms of a single parameter r\ — uji/lu c , 1/2 < rj < 1 
that describes the size of the photonic bandgap, giving: 
uj+ = u c , c + = c, 7 + = fe (l - rf), uj- = u> c {rf)/ (2T) - 1), 
c_ = c7]/^/2rj - 1, and 7_ = k (l - r\)j ^J2r] - 1. 



From the dispersion relation (Bl), the corresponding 
DOS, i. e., N(oj) = J d 3 kS(cj - uj{k)) is given by 



Aire 



2 [ fc o-\/0' 



N m (u) = < 



47rc: 



^/( W - W _)Vci-7i 

for < uj < u>i 

2 [fco+V^-^ + ) 2 /^- 7 ^] 2 (^-^ + ) 
V^-^+) 2 /c+- 2 -7i 

for uj u < uj < oo 



(B2) 



where a(0) = 1 and /3 M (0) = for all \x. 



For sufficiently large gaps (r\ < 0.9) and bare eigenfre- 
quencies lvq of the radiating dipole close to the upper 
band edge, it is an excellent approximation to ignore 
the lower branch of the photon dispersion relation, i. e., 
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for k < fco. The resulting DOS for this so-called three- 
dimensional isotropic, one-sided bandgap model is shown 
in Fig. 1 for a value of gap width parameter -q — 0.8 and 
the gap center frequency uj c a/2nc = 0.5. The square-root 
singularity at the band edge as well as the UV divergence 
N m (uj) oc uj 2 as u) — > oo are clearly visible. 
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Figure 1: The DOS for the three-dimensional, isotropic 
one-sided bandgap model of a PC. The parameters (see 
appendix B) are r\ = 0.8 and u c a/2-Kc — 0.5. 

Figure 2: The radiation dynamics resulting from the 
three-dimensional, isotropic one-sided bandgap model 
DOS as shown in Fig. 1 for various values of the bare 
dipole oscillator frequency loq relative to the upper pho- 
tonic bandedge lo u . The photonic bandedge is siutated 
at uj u a/2irc — 0.6 and the bare dipole oscillator frequen- 
cies are (a) woa/27rc = 0.58, (b) u>oa/27rc = 0.595, (c) 
ujQa/2-Kc = 0.599, (d) ujQa/2irc = 0.6, (c) u>oa/2irc — 
0.601, (f) w a/27rc = 0.605, and (g) w a/27rc = 0.62. 
Clearly visible are normal-mode oscillations, or vacuum 
Rabi oscillations, and the fractional localization of radia- 
tion near the photonic bandedge. The coupling strength 
has been chosen such that g(u)o) = 10~ 4 . For frequen- 
cies deep in the photonic bandgap (uJoa/2irc = 0.58) and 
deep in the photonic conduction band (u>oa/2Trc = 0.62) 
we observe negligible and exponential decay, respectively. 
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Figure 1: The DOS for the three-dimensional, isotropic one-sided bandgap 
model of a PC. The parameters (see appendix B) are r\ = 0.8 and Lo c a/2nc = 
0.5. 
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Figure 2: The radiation dynamics resulting from the three-dimensional, 
isotropic one-sided bandgap model DOS as shown in Fig. 1 for various 
values of the bare dipole oscillator frequency u)q relative to the upper pho- 
tonic bandedge co u . The photonic bandedge is siutated at LU u a/2nc = 0.6 and 
the bare dipole oscillator frequencies are (a) looo,/2ttc = 0.58, (b) looo,/2ttc = 
0.595, (c) LO a/27rc = 0.599, (d) co a/27rc = 0.6, (e) co a/27rc = 0.601, (f) 
ujQa/2TTc = 0.605, and (g) looo,/2ttc = 0.62. Clearly visible are normal- 
mode oscillations, or vacuum Rabi oscillations, and the fractional localiza- 
tion of radiation near the photonic bandedge. The coupling strength has 
been chosen such that g(uo) = 10~ 4 . For frequencies deep in the pho- 
tonic bandgap (iOoa/2Trc = 0.58) and deep in the photonic conduction band 
(u}Qa/2-Kc = 0.62) we observe negligible and exponential decay, respectively. 
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